In [1]:
import numpy as np
import matplotlib.pyplot as plt
from numba import njit,prange
import time
from matplotlib.cm import ScalarMappable
import matplotlib.animation as animation
from matplotlib.colors import LinearSegmentedColormap
import math
import warnings as kxe_kxe
kxe_kxe.filterwarnings('ignore')
In [2]:
@njit(nopython=True) # >10 min
def integ_u_t_nopar(x,y,V,D_x1,D_y1,I_u,x_ku,y_ku,sigma_U_x,sigma_U_y):
N,M = V.shape
out = 0
i_min = max(0, x - D_x1); i_max = min(M-1, x + D_x1)
j_min = max(0, y - D_y1); j_max = min(N-1, y + D_y1)
denominator_x = 2 * sigma_U_x**2
denominator_y = 2 * sigma_U_y**2
for i in range(i_min,i_max + 1):
for j in range(j_min,j_max + 1, nogil=True, reduction='+'):
out += I_u * np.exp( -1 * (i - x_ku)**2 / (denominator_x) - (j - y_ku)**2 / (denominator_y) ) * V[j,i]
return out
@njit(parallel=True) # 5:48 min
def integ_u_t_parallel_try(x,y,V,D_x1,D_y1,I_u,x_ku,y_ku,sigma_U_x,sigma_U_y):
N,M = V.shape
out = 0
i_min = max(0, x - D_x1); i_max = min(M-1, x + D_x1)
j_min = max(0, y - D_y1); j_max = min(N-1, y + D_y1)
denominator_x = 2 * sigma_U_x**2
denominator_y = 2 * sigma_U_y**2
for i in prange(i_min,i_max + 1):
for j in prange(j_min,j_max + 1):
out += I_u * np.exp( -1 * (i - x_ku)**2 / (denominator_x) - (j - y_ku)**2 / (denominator_y) ) * V[j,i]
return out
@njit(nopython=True, parallel=True) # 4:49 min
def integ_u_t_optimized(x, y, V, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y):
N, M = V.shape
out = 0.0
i_min = max(0, x - D_x1)
i_max = min(M - 1, x + D_x1)
j_min = max(0, y - D_y1)
j_max = min(N - 1, y + D_y1)
denom_x = 2.0 * sigma_U_x**2
denom_y = 2.0 * sigma_U_y**2
coeff_x = -1.0 / denom_x
coeff_y = -1.0 / denom_y
# Предвычисление экспонент для x
exp_x = [0.0] * (i_max - i_min + 1)
for idx_i, i in enumerate(range(i_min, i_max + 1)):
di = i - x_ku
exp_x[idx_i] = math.exp(di * di * coeff_x)
# Предвычисление экспонент для y
exp_y = [0.0] * (j_max - j_min + 1)
for idx_j, j in enumerate(range(j_min, j_max + 1)):
dj = j - y_ku
exp_y[idx_j] = math.exp(dj * dj * coeff_y)
# Параллельный цикл по j
total = 0.0
for idx_j in prange(len(exp_y)):
j_val = j_min + idx_j
ey = exp_y[idx_j]
local_sum = 0.0
for idx_i in range(len(exp_x)):
i_val = i_min + idx_i
ex = exp_x[idx_i]
local_sum += I_u * ex * ey * V[j_val, i_val]
total += local_sum
return total
integ_u_t = integ_u_t_optimized
In [ ]:
@njit
def T(x,y,V,U,b_v,d_v,c_u,integ):
return b_v * V[y,x] - d_v * V[y,x]**2 - c_u * U[y,x] * V[y,x] + integ
@njit(nopython=True)
def plots_single_try(h,U,V,
c_u,b_u,d_u,I_u,x_ku,y_ku,sigma_U_x,sigma_U_y,
c_v,b_v,d_v,I_v,x_kv,y_kv,sigma_V_x,sigma_V_y,
D_x1,D_y1,
D_x2,D_y2):
y_lim,x_lim = V.shape
U_new,V_new = np.zeros(U.shape),np.zeros(V.shape)
k0,k1,k2,k3 = np.zeros(V.shape),np.zeros(V.shape),np.zeros(V.shape),np.zeros(V.shape)
l0,l1,l2,l3 = np.zeros(U.shape),np.zeros(U.shape),np.zeros(U.shape),np.zeros(U.shape)
# k0
# l0
for x in range(x_lim):
for y in range(y_lim):
integ = integ_u_t(x,y,V,D_x1,D_y1,I_u,x_ku,y_ku,sigma_U_x,sigma_U_y)
k0[y,x] = T(x,y,V,U,b_v,d_v,c_u,integ)
integ = integ_u_t(x,y,U,D_x2,D_y2,I_v,x_kv,y_kv,sigma_V_x,sigma_V_y)
l0[y,x] = T(x,y,U,V,b_u,d_u,c_v,integ)
# k1
# l1
V1 = V + h/2 * k0
U1 = U + h/2 * l0
for x in range(x_lim):
for y in range(y_lim):
integ = integ_u_t(x,y,V1,D_x1,D_y1,I_u,x_ku,y_ku,sigma_U_x,sigma_U_y)
k1[y,x] = T(x,y,V1,U1,b_v,d_v,c_u,integ)
integ = integ_u_t(x,y,U1,D_x2,D_y2,I_v,x_kv,y_kv,sigma_V_x,sigma_V_y)
l1[y,x] = T(x,y,U1,V1,b_u,d_u,c_v,integ)
# k2 и l2
V2 = V + h/2 * k1
U2 = U + h/2 * l1
for x in range(x_lim):
for y in range(y_lim):
integ_v = integ_u_t(x, y, V2, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k2[y,x] = T(x, y, V2, U2, b_v, d_v, c_u, integ_v)
integ_u = integ_u_t(x, y, U2, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l2[y,x] = T(x, y, U2, V2, b_u, d_u, c_v, integ_u)
# k3 и l3
V3 = V + h * k2
U3 = U + h * l2
for x in range(x_lim):
for y in range(y_lim):
integ_v = integ_u_t(x, y, V3, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k3[y,x] = T(x, y, V3, U3, b_v, d_v, c_u, integ_v)
integ_u = integ_u_t(x, y, U3, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l3[y,x] = T(x, y, U3, V3, b_u, d_u, c_v, integ_u)
# Финальное обновление
V_new = V + (h/6) * (k0 + 2*k1 + 2*k2 + k3)
U_new = U + (h/6) * (l0 + 2*l1 + 2*l2 + l3)
return U_new, V_new
@njit(nopython=True, fastmath=True, parallel=True)
def plots_single_opt(h, U, V,
c_u, b_u, d_u, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y,
c_v, b_v, d_v, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y,
D_x1, D_y1,
D_x2, D_y2):
y_lim, x_lim = V.shape
U_new = np.zeros(U.shape)
V_new = np.zeros(V.shape)
k0 = np.zeros(V.shape)
k1 = np.zeros(V.shape)
k2 = np.zeros(V.shape)
k3 = np.zeros(V.shape)
l0 = np.zeros(U.shape)
l1 = np.zeros(U.shape)
l2 = np.zeros(U.shape)
l3 = np.zeros(U.shape)
h_over_2 = h / 2
h_over_6 = h / 6
# k0, l0
for x in prange(x_lim):
for y in range(y_lim):
integ = integ_u_t(x, y, V, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k0[y, x] = T(x, y, V, U, b_v, d_v, c_u, integ)
integ = integ_u_t(x, y, U, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l0[y, x] = T(x, y, U, V, b_u, d_u, c_v, integ)
# k1, l1
V1 = V + h_over_2 * k0
U1 = U + h_over_2 * l0
for x in prange(x_lim):
for y in range(y_lim):
integ = integ_u_t(x, y, V1, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k1[y, x] = T(x, y, V1, U1, b_v, d_v, c_u, integ)
integ = integ_u_t(x, y, U1, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l1[y, x] = T(x, y, U1, V1, b_u, d_u, c_v, integ)
# k2, l2
V2 = V + h_over_2 * k1
U2 = U + h_over_2 * l1
for x in prange(x_lim):
for y in range(y_lim):
integ_v = integ_u_t(x, y, V2, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k2[y, x] = T(x, y, V2, U2, b_v, d_v, c_u, integ_v)
integ_u = integ_u_t(x, y, U2, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l2[y, x] = T(x, y, U2, V2, b_u, d_u, c_v, integ_u)
# k3, l3
V3 = V + h * k2
U3 = U + h * l2
for x in prange(x_lim):
for y in range(y_lim):
integ_v = integ_u_t(x, y, V3, D_x1, D_y1, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y)
k3[y, x] = T(x, y, V3, U3, b_v, d_v, c_u, integ_v)
integ_u = integ_u_t(x, y, U3, D_x2, D_y2, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y)
l3[y, x] = T(x, y, U3, V3, b_u, d_u, c_v, integ_u)
# Финальное обновление
V_new[:] = V + h_over_6 * (k0 + 2*k1 + 2*k2 + k3)
U_new[:] = U + h_over_6 * (l0 + 2*l1 + 2*l2 + l3)
return U_new, V_new
plots_single = plots_single_opt
Для примера пары конкурирующих популяций можно рассмотреть сардину-иваси (Sardinops sagax) и японского анчоуса (Engraulis japonicus), которые конкурируют за общий ресурс — планктон. При высокой численности сардин-иваси снижается популяция японского анчоуса из-за ограниченного доступа к пищевым ресурсам.
Обоснование:¶
- Межвидовая конкуренция: Оба вида питаются планктоном, что делает их прямыми конкурентами за пищевой ресурс.
- Эксплуатационная конкуренция: Сардины-иваси эффективнее потребляют планктон, что приводит к снижению его доступности для анчоусов .
- Влияние на численность: При увеличении численности сардин-иваси наблюдается уменьшение численности японского анчоуса, что подтверждает негативное влияние одной популяции на другую.
Этот пример иллюстрирует типичную косвенную конкуренцию (эксплуатационную), когда организмы конкурируют за общий ресурс, не взаимодействуя напрямую
In [48]:
from PIL import Image
img1 = Image.open("data\Sardinops_melanostictus.jpg")
img_array1 = np.array(img1, dtype=float) / 255.0
img2 = Image.open("data\Engraulis_japonicus_01.jpg")
img_array2 = np.array(img2, dtype=float) / 255.0
fig, axes = plt.subplots(1, 2, figsize=(10, 10))
axes[0].imshow(img_array2)
axes[0].axis('off')
axes[0].set_title('Engraulis japonicus')
axes[1].imshow(img_array1)
axes[1].axis('off')
axes[1].set_title('Sardinops sagax')
plt.tight_layout()
plt.show()
In [62]:
def make_gif(h,
c_u, b_u, d_u, I_u, sigma_U_x, sigma_U_y,
c_v, b_v, d_v, I_v, sigma_V_x, sigma_V_y,
D_x1, D_y1,
D_x2, D_y2,
num_frames=100, fps=10, filename='population_dynamics.gif', label1='Популяция_1', label2='Популяция_2', gridsize=100):
start_time = time.time()
x_ku = y_ku = 0.8*gridsize
x_kv = y_kv = 0.2*gridsize
U0 = np.zeros((gridsize, gridsize))
V0 = np.zeros((gridsize, gridsize))
# Начальные популяции в углах
U0[ gridsize*1//8: gridsize*3//8, gridsize*1//8: gridsize*3//8] = 1.0 # Левый верхний угол
V0[-gridsize*3//8:-gridsize*1//8, -gridsize*3//8:-gridsize*1//8] = 1.0 # Правый нижний угол
U = U0.copy()
V = V0.copy()
# Фиксированные цветовые карты
colors_v = [(0,0,0,0), (1,1,1,0.3), (1,0,0,1)] # Градиент от белого до красного
colors_u = [(0,0,0,0), (1,1,1,0.3), (0,0,1,1)] # Градиент от белого до синего
cmap_v = LinearSegmentedColormap.from_list(label1, colors_v)
cmap_u = LinearSegmentedColormap.from_list(label2, colors_u)
# Фиксированные пределы нормализации (если известны)
vmax_v = max(1e-3, np.max(V0))
vmax_u = max(1e-3, np.max(U0))
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_title("Динамика конкурирующих популяций")
ax.set_xticks([])
ax.set_yticks([])
img_v = ax.imshow(V0, cmap=cmap_v, vmin=0, vmax=vmax_v, alpha=0.7, interpolation='bilinear')
img_u = ax.imshow(U0, cmap=cmap_u, vmin=0, vmax=vmax_u, alpha=0.7, interpolation='bilinear')
cbar_v = fig.colorbar(ScalarMappable(cmap=cmap_v), ax=ax, label=label1)
cbar_u = fig.colorbar(ScalarMappable(cmap=cmap_u), ax=ax, label=label2)
def update(frame):
nonlocal U, V
U, V = plots_single(h, U, V,
c_u, b_u, d_u, I_u, x_ku, y_ku, sigma_U_x, sigma_U_y,
c_v, b_v, d_v, I_v, x_kv, y_kv, sigma_V_x, sigma_V_y,
D_x1, D_y1,
D_x2, D_y2)
U[U < 0] = 0
V[V < 0] = 0
img_v.set_array(V)
img_u.set_array(U)
ax.set_title(f"Шаг {frame}, t={frame*h:.2f}")
return img_v, img_u
ani = animation.FuncAnimation(fig, update, frames=num_frames, interval=1000//fps, blit=True)
ani.save(filename, writer='pillow', fps=fps, dpi=100)
plt.close()
print(f"Анимация сохранена как {filename}")
print(f'На выполнение расчетов понадобилось {(time.time() - start_time):.4f} секунд')
return ani
# ani = make_gif(h=0.05,
# c_u=0.05, b_u=0.7, d_u=0.15, I_u=0.003, sigma_U_x=20, sigma_U_y=20,
# c_v=0.1, b_v=0.8, d_v=0.1, I_v=0.003, sigma_V_x=20, sigma_V_y=20,
# D_x1=30, D_y1=30,
# D_x2=30, D_y2=30,
# num_frames=2000, fps=25, filename='populations.gif',gridsize=800, label1='Engraulis japonicus',label2='Sardinops sagax')
ani = make_gif(h=0.05,
c_u=0.05, b_u=0.7, d_u=0.15, I_u=0.003, sigma_U_x=20, sigma_U_y=20,
c_v=0.1, b_v=0.8, d_v=0.1, I_v=0.003, sigma_V_x=20, sigma_V_y=20,
D_x1=30, D_y1=30,
D_x2=30, D_y2=30,
num_frames=100, fps=10, filename='populations.gif',gridsize=100, label1='Engraulis japonicus',label2='Sardinops sagax')
Анимация сохранена как populations.gif На выполнение расчетов понадобилось 26.2042 секунд
In [63]:
from IPython.display import HTML
HTML(ani.to_jshtml())
Out[63]: